A theoretical approach to zero-reflection toroidal curved metasurfaces

In this paper, we investigate the electromagnetic response of metasurfaces due to excitation of the toroidal moment. A toroidal curved metasurface analyzad using a novel theoretical solution based on the Fourier analysis to evaluate the localized fields. Analyzing localized near-field interactions are crucial in investigating the excited trapped modes and enables us to optimize the reflection properties of the proposed metasurface. Optimization is accomplished using graphene layer and resulted a hybrid dielectric-graphene structure with near-zero reflection properties.

A theoretical approach to zero-reflection toroidal curved metasurfaces Hosein Allahverdizade 1 , Sina Aghdasinia 1 , Hemn Younesiraad 1,2 & Mohammad Bemani 1* In this paper, we investigate the electromagnetic response of metasurfaces due to excitation of the toroidal moment. A toroidal curved metasurface analyzad using a novel theoretical solution based on the Fourier analysis to evaluate the localized fields. Analyzing localized near-field interactions are crucial in investigating the excited trapped modes and enables us to optimize the reflection properties of the proposed metasurface. Optimization is accomplished using graphene layer and resulted a hybrid dielectric-graphene structure with near-zero reflection properties.
Approaches for calculating scattered electromagnetic fields of metasurfaces have been proposed in the past, such as using periodic Green's function (Floquet-Bloch) or the GSTC method 1 , but the fundamental issue with applying these methods is that the structures must have a minimum of periodicity in order to be employed and they are unsuitable for calculating near-field responses.
A zero-reflection structure is required to achieve transparency. This characteristic is present in the wings of the remarkable species of butterfly known as Greta oto Fig. 1a 2,3 . According to a recent study 3 , Greta Otto's wings have asymmetrically dispersed nanostructures in the transparent areas Fig. 1b. It's interesting that this aspect of the butterfly's wings is what makes it transparent from different angles 4 . This observation motivates us to investigate analytical approaches of non-periodic nano-structures.
The electromagnetic response of metasurfaces can be related to traditional electric and magnetic dipoles or their complex combinations known as multipoles 5 . A toroidal dipole is the third member of localized electromagnetic excitations which is created by a current circulating on the surface of a torus and has been observed in solid-state physics 6 . The toroidal response of metasurfaces was experimentally observed in microwave regime 7 and then theoretically scaled to the THz regime 8 . The performances of these metasurfaces are usually limited by radiative and nonradiative losses where nonradiative losses can by reduced by employing materials of low loss such as dielectric 9,10 likewise the solution employed in nontoroidal metasurfaces 11 .
Radiative losses can be regulated by designing the asymmetric all-dielectric unit cells to excite a high-quality factor resonant response known as Fano resonance in which the radiative damping can be efficiently suppressed by trapped mode and leads to the reduction of radiative losses 12 . Trapped modes that excited in such asymmetrical structures are included in the the concept of bound states in the continuum (BIC) 13 . A strong link between the toroidal dipole resonance and the BIC was defined in the context of all-dielectric metasurfaces 14 . To construct a trapped mode supporting metasurface, we utilize a set of identical subwavelength unit cells with two asymmetric high-refractive index Silicon particles made in the form of modified nanodisks.
Trapped modes lead to a strong near-field enhancement where the located electromagnetic fields becomes notable. We utilize a novel theoretical approach to analytically evaluate localized electromagnetic fields near the all-dielectric unit cells of trapped mode supporting metasurface.
This evaluation provides the capability of manipulating the reflection coefficient for the surface. Our investigations should support every possible surface shape in the practical applications of such metasurfaces as sensing 15,16 . To accomplish a comprehensive solution, we consider a cylindrical surface where any arbitrary surface can be fitted on a cylindrical surface with a defined cross-section.
Since all-dielectric curved metasurfaces are composed of multi subwavelength cells, their optical properties are related to resonant features of each constitutive unit cell resonator and its mutual interaction with the other cells, trapped mode excitation does not merely satisfy the high-Q and zero-reflection properties which necessitates computing the localized electromagnetic fields to optimize such curved multi-cell metasurfaces.
The sampled curved metasurface was selected as ρ surface in Fig. 2 with complete ϕ angle which enables us to analyze the localized electromagnetic fields with Fourier series analysis in 2π interval which is discreted to 20 individual values to easily manipulate the sections to reduce the overall reflectivity.
Graphene is selected to optimize the obtained reflectivity from evaluated electromagnetic fields. Optimization can be implemented by varying the graphene chemical potential when the graphene layer is integrated with the proposed dielectric structure.

Results and methods
Theoretical formulation. The toroidal moment T emerges from the multipole expansion of an arbitrary localized current distribution 6 . T is defined in terms of current density distribution j(r) , as follows where c indicates the speed of light in free space. Equation (1) is satisfied by the following form of current density distribution By decomposing current vector into longitudinal j and transversal j ⊥ parts, Eq. (1) can be recast into a more convenient form. j does not contribute to the toroidal moment T 6 but j ⊥ can be written as the curl of the magnetization density j ⊥ = c∇ × µ(r) and gives the toroidal moment in terms of magnetization density 17 (1) T = 1 10c r(r.j) − 2r 2 j d 3 r (2) j(r) = c∇ × ∇ × δ(r)T.  www.nature.com/scientificreports/ Equation (3) can be used to evaluate the toroidal moment for a finite system with known magnetization density. Representing the magnetization density by a distribution of localized magnetic moments m α at sites r α leads to the following equation for localized magnetization density µ loc hence the toroidal moment takes the following form simplified to localized magnetic moments Figure 3 shows a toroidal construction made up of two cylinders of various sizes. Based on Eq. (5), this structure clearly illustrates how to distribute magnetic moments in two opposite directions and so form a toroidal moment. Because of the simplicity of modeling and the ability of application to different forms figure, we will utilize this structure in our investigations 18 .
Respectively the potential caused by the toroidal moment A T takes the following form due to the similarity of T and the electrical moment P 19 The localized magnetic field at arbitrary sites, H loc is equal to the field caused by an external source, H ext and integral of the indicated potential caused by toroidal moment however the toroidal moment is itself affected by external source. We know that e jk|r−r ′ | 4π|r−r ′ | term in Eq. (6) is the free-space Green function in the spherical coordinates 20 which can be replaced with its equivalent in the cylindrical coordinates www.nature.com/scientificreports/ where H 2 0 (ρ) is Hankel function of the second kind, zeroth order and by assuming H ext in the cylindrical coordinates, the total localized magnetic field can be defined in cylindrical coordinates which enables us to analyze surfaces like ρ surface in Fig. 1 where ρ is dependent to ϕ angle.
The toroidal moment can be defined on each unit cell of ρ surface produced by non-directional magnetic moments of two nanodisks with d spacing due to 2π period of ρ surface, H loc can be expanded in the Fourier series with H loc, coefficients as follows evaluating the toroidal potential caused by this expanded localized H-field, leads to calculating overall H loc in Eq. (7) by a novel approach in Methods. Another double Fourier series with H s,x l,l ′ , H s,y l,l ′ and H s,z l,l ′ coefficients appears in calculated H loc as follows where α m produced as we used a linear relationship between the magnetic field and the resulting magnetization which can be obtained in structures similar to a cylinder with L length and 2r 1 diameter as 21 where n s = 3.88 is the index of refraction of the silicon nanodisk medium 22 . However our medium has been modified by embedding two extra 30 • inclined nonodisks with respect to main nanodisk as Fig. 1 to make the response stable for oblique incidents which is compared for 2 different external incidents in Fig. 4.

Methods.
This section describes how to compute a localized magnetic field using the previously mentioned potential formula (8). Since the localized magnetic field has a Fourier series form because of the 2π periodicity in our structure, we can use some basic Fourier series properties to solve the integral equation. Replacing Hankel part of G by double Fourier series with a l,l ′ coefficient and substituting Eq. (9) in (8), overall localized H-field equals where ( â ϕ ′ × H loc ) can be calculated as and results with respect to Fourier properties  (13) results Based on the Fourier series properties, multiplication of two periodic terms in (16) results in an integrated double Fourier series in which the coefficients are the convulsions of a l,l ′ with b l ′s summation on l is independant of ( ∇ ′ × ) operator and d ′ integration. ( ∇ ′ × ) operates on e jl ′ ϕ ′â x , e jl ′ ϕ ′â y and e jl ′ ϕ ′â z which can for example be calculated as ∇ ′ × (e jl ′ ϕ ′â z ) = jl ′ ρ ′ e jl ′ ϕ ′â ρ ′ where â ρ ′ will be replaced with equivalent Cartesian unit vectors. Eventually (17) takes the following form Due to the mentioned singularity in evaluating localized field, angular region of singularity �ϕ = �P/nd where nd is the approximate radial distance of unit cell and n is in the order of 1000, will be subtracted from the integration region so Another Fourier series with d l ′ coefficient assumed as follows multiplying the mentioned Fourier series by the other items in (18) summarized to a single Fourier series where For the convolution operations we pursue the following approach 23 which reduces all the convolutions to the multiplication of coefficients as follows Now according to the Fourier series properties, following integral of ϕ ′ in 2π interval reduces to a constant item of the corresponding series where u ′ = 0 and we have (15)  www.nature.com/scientificreports/ Analysis verification. Excitation and surface reflectivity. The theory presented in previous section can be used to calculate the localized field on surfaces like ρ in the presence of appropriate external field which leads to obtaining reflection coefficient on the surface. Line source I 0 at (ρ 0 , ϕ 0 ) with ω = 2π × 2.67 THz considered as an external source 24 proved that this source produces a TM z wave in ρ ≪ ρ o with the following E-field which is in the form of (8) and can be expanded in the Fourier series, Hankel part of E ext z replaced by a Fourier series with f l coefficient and the corresponding transverse H-field components can be calculated as 24 which results by substituting cos ϕ ρ by a series with g ′ l coefficient and sin ϕ ρ by a series with g ′′ l , we have the x component of external H-field as The normalized x component of external H-field is obtained as another Fourier series with H ext,x l coefficient which can be evaluated at 20 discrete values of ϕ on ρ surface where the surface was selected as an elliptical cylinder with cross-section equation ρ = 2/ 3 cos 2 ϕ + 1 . The evaluated Fourier series at each determined discrete value using MATLAB software demonstrated in Fig. 5.
The overall x component of localized H-field in the presence of proposed external source is obtained using (11) as and demonstrated in Fig. 6 at the same 20 discrete values of ϕ angle. Same procedure could be used to evaluate another series for H loc y respectively, but the selected external source causes H loc y ≪ H loc x at ρ ≪ ρ o and makes the y component of reflection coefficient less important.   Graphene strip lines. Graphene Strip Lines are added as patch arrays to reduce the overall reflectivity of proposed metasurface which is examined for a unit cell in Fig. 9. It is clearly obvious that strip lines reduces overall reflectivity of our structure and eliminates the real part of impedance. Focusing on eliminating the imaginary part of impedance is beneficial in achieving the zero-reflection properties for the surface in the following. The impedance of this layer for TM-polarized incident field is determined as 26 where k 0 is the wave number in free space. Host medium incident wave number k eff and wave impedance of that η eff are obtained from relative effective permittivity ε eff = ε r +1 2 where ε r in the desired graphene medium is defined as 27 where the scalar surface conductivity of layer is σ g = σ intra + σ inter in which the former term is due to the intraband contributions and the latter term is due to the interband contributions. In mid-infrared frequencies where ω is below 2µ c / , the interband term is negligible and the intraband term is dominant and equals 28 where Ŵ is a phenomenological scattering rate, k B is Boltzmann's constant, is the reduced Planck constant and T is the temperature. θ in Eq. (32) is the angle of incident vector with ρ surface which can be written in the terms of ϕ angle as follows and the grid parameter α in Eq. (32) is defined as 25 where the strip width and their spacing are selected as W = 2.58 µm , D = 5.16 µm for 2.67 THz frequency. It is crucial to note that every optimization involves making a trade-off. For instance, the use of strip lines can reduce reflectivity, but it also induces capacitive and inductive effects on the structure 29 , leading to undesired outcomes like frequency shifts in the main frequency.
Graphene ribbons. In this section, we employ a knock-out technique to achieve our ultimate goal of zeroreflection. By applying different potentials to graphene ribbons, we can attain the desired surface impedance, which can be manipulated to achieve zero-reflection properties.
Graphene ribbons with t g = 10 nm thickness and distinct values of chemical potential µ c are added to each interval of ϕ angle, n. Z g is caused by this layer Figure 9. A unit cell reflectivity properties with applying graphene strip lines and without utilizing it. www.nature.com/scientificreports/ The overall total impedance of designed 3-layer hybrid structure is 30 the total reflectivity of the structure R t , obtained from the overall total impedance, is dependant to µ c of the selected ribbons in each interval of ϕ angle, n. The selection procedure will done by observing the imaginary part of R t in the presence of 11 different quantities of µ c in [0, 1] interval of it, denoted by n ′ . The imaginary Figure 10. The imaginary part of R t in 4 different intervals of ϕ angle. www.nature.com/scientificreports/ part of conductivity plays an important role in the propagation of surface waves guided by the graphene sheet 31 so the optimizing procedure is restricted to cancelling the imaginary part of R t . Figure 10 shows the imaginary part of R t in 4 different intervals of ϕ angle, n for 11 possible selections of µ c , n ′ . All the optimized µ c selections are shown in Table 1.
Finally a reasonable segment of the optimized 3-layer hybrid structure involving five intervals of ϕ angle from n = 8 to n = 12 as Fig. 11 is simulated in CST Microwave Studio to validate the accomplished optimization procedure. In non-optimized case all the graphene ribbons are simulated with µ c = 0 . Excited modes are obtained due to the Eigenmode Solver for both optimized and non-optimized structures and the reflectivity |R| in corresponding wavelengths are compared in each mode as Fig. 12.
By applying different chemical potentials to graphene ribbons and achieving varying conductivities, excitation modes may be formed with slight shifts in their frequencies. For instance, the second mode of the non-optimized structure appears at 2.64 THz with an equivalent reflectivity of approximately 0.4, while the second mode of the optimized structure appears at 2.67 THz with an equivalent reflectivity near to zero. The S-parameter and the reflectivity |R| are showed in details as Fig. 13.

Discussion
In summary, we introduced a novel approach to determine the electromagnetic response of dielectric-graphene hybrid curved metasurfaces due to excitation of the toroidal moment. Analyzed metasurface with asymmetric unit cells designed for exciting the well-known trapped mode whereas they show a robust response for oblique incidents. Trapped mode excitation and the concept of bound states in the continuum (BIC) do not merely provide the high quality properties for proposed curved structure and accentuates our introduced novel approach for determining localized electromagnetic fields. Consequently, any practical surface with random and arbitrary curve can be optimized by utilizing the near-field responses. To figure out the functionality of the procedure, a reasonable segment of the designed Silicon metasurface is optimized using graphene ribbons and strip lines which results near-zero reflectivity properties for the structure. To bias graphene sheets empirically, one can apply an external voltage to the sheets using a voltage source or bias tee. Graphene can be indirectly stimulated by applying thin gold electrodes on the side borders of the structure and mounting the graphene layers onto them. By applying stimulation to the metal electrodes, the graphene layers can be activated or manipulated to achieve the desired conductivity or impedance properties. It is important to note that the exact procedure for mounting and stimulating graphene layers may depend on the specific experimental setup and conditions, and may require careful calibration and optimization to achieve the desired results 32 . In a more advanced approach, tension-sensitive resistors can be utilized in a way that allows the change of curvature to result in a variable  www.nature.com/scientificreports/ voltage. This can enable the structure to self-adjust, as the resistors can dynamically respond to changes in the curvature or deformation of the structure.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.